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ABSTRACT 

The one-dimensional, ordinary differential equation (ODE) that satisfies the midplane gravitational potential of truncated, flat power- 
law disks is extended to the whole physical space. It is shown that thickness effects (i.e. non-flatness) can be easily accounted for by 
implementing an appropriate "softening length" A. The solution of this "softened ODE" has the following properties: i) it is regular 
at the edges (finite radial accelerations), ii) it possesses the correct long-range properties, iii) it matches the Newtonian potential 
of a geometrically thin disk veiy well, and iv) it tends continuously to the flat disk solution in the limit /I — » 0. As illustrated by 
many examples, the ODE, subject to exact Dirichlet conditions, can be solved numerically with efficiency for any given colatitude at 
second-order from center to infinity using radial mapping. This approach is therefore particularly well-suited to generating grids of 
gravitational forces in order to study particles moving under the field of a gravitating disk as found in various contexts (active nuclei, 
stellar systems, young stellar objects). Extension to non-power-law surface density profiles is straightforward through superposition. 
Grids can be produced upon request. 

Key words. Accretion, accretion disks — Gravitation — Methods: analytical — Methods: numerical 
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1. Introduction 

Flattened astrophysical objects like disks are known to gen- 
rally produce gravitational fields weaker than spherical bodies 
of comparable mass. However, the gravity of low-mass disks 
evolving on long time scales may play a significant role in 
their own dynamics and environment (Goldreich & Tr emainel 
[19781: ISubr&Karasll2005h . Gravitational forces from disks are 
not easily accessible by numerical computation, and this do- 
main still represents an interesting challenge. For various rea- 
sons (misknowledge of boundary conditions, sensitivity and in- 
accuracy of solutions, relevant physical scales, kernel singulari- 
ties, computational cost, etc.), neither the Poisson equation nor 
Newton's integral law offers a simple and straightforward tool, 
and each must be handled with some caution. Truncated expan- 
sio ns of solutions generally suffer from inaccuracy and instabil- 
ity jClem ent 1974; Hachisu 1986). Softened Gravity for con- 
tinuous systems may help in some circumstances, but the influ- 
ence of the softening length — a free-parameter, classically — 
is spurious, and it fundamentally destroys the Newtonian charac- 
ter of the gra vitational interaction (iHocknev & Eastwood 1988; 
lAdams et alJ[l989) . Each disk configuration (symmetry, edges, 
mass profile, shape, etc.) must therefore be investigated individ- 
ually for a given application. 

Geometrically thin disks probably constitute the main class 
of astrophysical disks. These exhibit various shapes, density pro- 
files and sizes. For those orbiting a massive central object (star or 
black hole), a self-similar behavior may develop secularly, lead- 
ing to a mass density profile varying close to a power law of the 
radius. Such a profile is widely supported by theory an d it is a 
typical initial ingr edient of numerical simulations (e.g. iPringld 
ll98lHEdgaiil2007h . Even in quasi-Keplerian rotation, thin disks 
can be influenced by their own gravity. Hure & Hersant (2007) 
(hereafter Paper I) have shown that the midplane gravitational 
potential of flat, power-law disks obey an ordinary differential 



equation (ODE) accounting for edges which are usually ignored 
(Bisnovatvi-Kogan 1975; Goodman & Evans 1999). Analytical 
solutions i n the form o f very rapidly converging serie s have been 
report ed in lHure et al.l (12008) . At the same time, Hu re & PierensI 
(l2009h have shown that the model of "Softened Gravity" offers a 
good framework for determining the Newtonian potential of thin 
disks (whatever the mass density profile), provided the "soften- 
ing length" takes a very specific form, locally. In the present pa- 
per, we s how that the ODE for t he gravitational potential de- 
scribed bv lHure & Hersani (l2007h can be extended to the whole 
physical space, and we combine this result with an appropriate 
softening length to describe the potential of disks of non vanish- 
ing thickness. 

The paper is organized as follows. We recall the basic con- 
figuration and useful formulae of potential in flat, power-law 
disks in Sect. [31 The derivation of the generalized ODE, non- 
dimensionalization and asymptotic behavior of solutions, are 
found in Sect. [3] The numerical solutions are given in Sect. [4] 
Thickness effects, including the introduction of the softening 
length, are discussed in Sect [5] The last section is devoted to 
a conclusion. 



2. Theoretical grounds and notation 

Following iHure & Hersanti (l2007l) (hereafter Paper I), we con- 
sider a flat axisymmetrical disk with inner edge ai„ > 0, outer 
edge flout > flin (see Fig. [Hi, and a power-law surface density of 
the form 



S(a) = 5:o - 

\flO 



(1) 



where a is the cylindrical radius, ao some reference radius, and 
So = S(flo) the corresponding surface density. This profile can 
serve as a basis for defining more complex mass distributions. 
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by mixing power laws with different indices (positive and nega- 
tive). For such a disk, the gravitational potential in space is given 
exactly by the expression (Durand 1953.) 



(A(r) 



where 



-2G 



I.(a)kK(k)da, 



Kik) 



Jo 



is the complete elliptic integral of the first kind, and 
2^/aR 



(2) 



(3) 



(4) 



is the modulus (0 < ^ < \), R and Z are the cylindrical coordi- 
nates, and r is the spherical radius (i.e. - + T?). Known 
properties about this configuration are the followings. The inte- 
gral in Eq. |2] has a diverging kernel as soon as the modulus k 
reaches unity. This occurs everywhere inside the disk. Standard 
quadrature schemes fail to give acc urate potential values u nless 
a specific treatment is considered (iHure &Pierensll2005l) . The 
potential is not a power law of the radius, because of edges. A 
closed form f or Eq. [2] exists only in the case of infinitely ex- 
tended disks (lBisnovatvi-Koganll975tlGoodman & Evansj 19991) 
and for finite size disks wi th constant surface density (iDurandl 
[1951 11 ass & Blitzeilll983h . The potential is finite everywhere, 
except when - and s + \ < 0, but has an infinite radial 
gradient in the r nidplane at edge cro ssing, as a consequence of 
flatness (Dur andl[T953t iMestell 119631) . Finally, the disk mass is 
finite except when - and i H- 2 < 0, or when aout °° and 
s -H 2 > 0. 




Fig. 1. Configuration for the flat, finite size disk with radius at 
the edges Oin and Oout; r and Q are spherical coordinates, R and Z 
are cylindrical coordinates. 



It has been established in Paper I that the potential \p given by 
Eq. |2]obeys, in the disk rnidplane (i.e. for Z = 0), a differential 
equation of the form 

— -(l-i-i)--A = 0, (5) 
dR R 

where A is fully analytical. Sol utions in th e form of rapidly con- 
verging series were reported in lHure et al.l(l2008l) . 



3. Generalized ODE 

The starting point of the generalization of Eq.|5]is the equation of 
the line containing the origin and making an angle Q G [0, «■] (the 
colatitude) with the z-axis (see Fig.[T]). This equation is simply 

Z = /?cotan0. (6) 



By setting u-ajR'm Eq.|4l and using Eq.|6] we get the relation 



- 2m I — - 1 1 + 1 + cotan^e = 0. 



(7) 



We then see that, once Q is specified, m is a function of the mod- 
ulus k only. This property was the condition for the existence of 
Eq. |5] It still holds here. We can therefore proceed as in Paper 
I : the partial derivative of the potential with respect to the ra- 
dius can be expressed as a function of the potential itself. Here, 
we consider the derivative with respect to the spherical radius r 
instead of the cylindrical radius R (although this choice is not re- 
ally important). The full demonstration is given in Appendix lAl 
Strictly speaking, we obtain a partial differential equation (PDE) 
that becomes, for a given colatitude 9, an ordinary differential 
equation (ODE). It writes as 



^-(1 + .)!^ 

dr r 



A = 0, 
where A is defined by 



A = 2G£() — M, 



"^^O -(j+l) 



"0 



(8) 



(9) 



with Mout = aout/R, and Mi,, = ain/R. This generalizes to the whole 
space the differential equation reported in Paper I, which is valid 
only in the disk plane (i.e. the case = |). The second member A 
is an analytical function of a single spatial coordinate r. It now 
depends on the colatitude 6, and on the four parameters s, aq, 
flin, and flout- The main difference comes from the Z-dependent 
modulus k. 

We can make the ODE scale-free by setting r - r/ao, 
and if/ = lA/iAo where i/'o is a constant. Then, we get the non- 
dimensional ODE: 



dr r 
where 



(10) 



(11) 



Potential values are known at the two edges of th e disk in 
the form of rapidly converging series (iHure et al ]l2008l) . and can 
therefore be used for adimensioning. A more convenient value 
for i/'o is probably the potential at the origin which is exact. From 
Eq. |2] we have 

iJ/iO) — — 27rGS()floMQ**^'^ Jj""' mWm 

= -2nG'Loaou^^'^^^u'^'^IXs+i = ^c, 
where we have set (definition different than in Paper I) 
1 - A" 



(12) 



Xn 



(13) 



and A - flin /flout is the axis ratio of the disk. Thus, with i/'o = iff^, 
the inhomogeneous term in the non-dimensional ODE is 



S = ^— r f V«^^outK(feout) - A'+' ^/u~kinK(ki„)\ . 



(14) 



When flin = (no inner edge; the disk extends down to the cen- 
ter), the last term vanishes. When flout — » this is the first 
term. The asymptotic properties of the solution ifr are discussed 
in Appendix IB] 



' The choice oq = " in (or Oq = aaai) is possible but not appropriate in 
the particular case where a[„ = (respect, flout — > °°)- 
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Fig. 2. Potential versus ( due to a flat disk with A - 0.1, 
flout/flo - 1 and power index s - -1.5, for several colatitudes 
9. Values are normalized to the central value i^c- The grid is such 
that A - n/4 and - 1000. Infinity stands at ^ = 2, and the 
edges are located at ^ » 0.1269 and 1 (for = |). Although not 
visible on the figure, the potential is not derivable at both edges 
(case with = see Sect.|2]i. 



Actually, a centered second-order space discretization based on 
+ 1 mesh points yields the system: 



<Ao = 1, 



iA„-i + 2f /„(1 + i)(A„ - + 2ff„Q„ = 0, 



« e {1,. 



.,A^- 1) 



(19) 



lAiv =0. 



with ^0 = 0, ^jv = ^, and ^„ = nd^. For s = -1, it can be 
advantageous to rewrite the ODE in a slightly diff^erent form, for 
instance by using the auxiliary function M = ijfr, which satisfies 
the following ODE: 



dM M 1 
{2 + s)f 

di: ^ C cos2(A0 



e = o. 



(20) 



with the boundary conditions M(0) - and m(^) = 

To illustrate this numerical part, we show in figure |2]the so- 
lution {(^,-, i/t,)) obtained at a few colatitudes with - 1000 and 
A - J (see note|2]i and for the following parameters: A = 0.1, 
flout/flo = 1, and s - -1.5 (due to adimensioning, it is not nec- 
essary to specify Eo). Figures[3]to|5]correspond to power indices 
s = {-1,-0.5,0). 



4. Numerical solution in space with mapping 

There are many methods for solving a first-order ODE numer- 
ically. This needs one boundary condition for an initial value 
problem (IVP), but two conditions for a two boundary value 
problem (TBVP). The potential is known at four places: at the 
origin, at the two edges, and at infinity. The origin is particularly 
well suited to finding the potential in space in a given direction r, 
as the corresponding boundary condition is the same regardless 
of the colatitude. Values at the disk edges are useful mainly for 
determining the potential inside the disk or in its neighborhood. 
For some applications, potential values are required at large dis- 
tances from the disk (e.g. 'Subr et al ] l2004tlSubr & Karasl l2005h . 
In this case, it can be interesting to use the trivial boundary con- 
dition at infinity ((/»—> as r —> oo). This is easily performed by 
mapping the whole space. For instance, the transformatioifl 



^(r) = — atanr. 



(15) 



where A is a positive constant, maps the range re [0, oo[ into the 
compact domain ^ e [0, ^]. The "mapped ODE" reads 



d\i/ if' Q 

-^-(l+.)/^-/f =0, 
d( i C 

where we have set 



- = sinc(2A0, 
and 



(16) 



(17) 



(18) 



It is then easy to discretize Eq. [16] on a grid, regular in ^, 
and to solve the associated linear system by standard techniques. 



^ Since the tangent function is essentially linear with the argument in 
the range [0, 7r/4], the choice A = n/4 leaves the space inside the shell 
r < flo almost undistorted by this change of coordinate. 




Fig. 3. Same legend as for Fig. |2]but for s 
used the auxiliary function M and Eql20l 



The time required to solve the above system is linear with 
A^. A full grid is then generated in a time that is proportional to 
NxM, where M is the number of mesh points in the 0-direction. 
With a classical laptop, this takes about 15s to generate a grid 
with A^^ = 10^ mesh points (including the update of the inho- 
mogeneous term Q). This is much faster than what can be ob- 
tained from the multipole expansion, by a factor proportional to 
the number of Legendre polynomials to be considered, which 
is gene rally very large — typically, several hundreds inside the 
source (IClementlll974 . 

Accuracy can be properly checked in the ho mogeneous case 
since the potential i s known in a closed form (lDurandlll953t 
lEass & Blitzei]fT98l . Figure |6] displays the error relative to the 
exact value. The radial grid and the disk parameters are the same 
as for Fig. |5](we used M = 100 for the sampling in 6). We see 
that the agreement is very good in the whole space, less than 
10"^ on average for this radial resolution, and it is rather uni- 
form in the whole space (however, with a certain deterioration 
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Fig. 4. Same legend as for Fig. |2]but for s 
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Fig. 5. Same legend as for Fig.|2]but for 5 = corresponding to 
an homogeneous disk. 



in the midplane due to edges; see Sect. |2]). The accuracy only 
depends on the discretization of the ODE, that is, only on (not 
on M). Increasing obviously improves potential values. 



5. Thickness effects and softening length 

The flat disk hypothesis is well-suited to many studies, both the- 
oretical and numerical ones. Its physical realism is, however, 
limited. For instance, gas or particles present at each edge would 
feel an infinite acceleration there (see again Sect. |2). Such sin- 
gularity would disappear by considering an extra dimension (or 
with a vanishing density profile at the disk edges). Thickness 
effects (matter present off" the midplane) can therefore be im- 
portant, not only for thermodynamic reasons, but also from a 
dynamical point of view. We have not yet investigated the exis- 
tence of an ODE corresponding to geometrically thin (i.e. non- 
flat) disks. However, we can reproduce the Newtonian potential 
due to a thin disk with a one dimensional differential equation 
of the form of E g. [8] by a "softened" potentiafl In this context, 
iHure & PierensI ((2009) showed that the midplane potential of a 
vertically stratified disk can be approximated by an equation re- 



^ For the concept of Soft ened Gra vity in point mass systems, see 
iHocknev & Eastwood ( 119881) , and see lAdams et d] ( fT989l) for the case 
of gas disks. 




-2 -1.5 



-0.5 



0.5 



1.5 



Fig. 6. Logarithm of the eiTor relative to exact values for the ho- 
mogeneous case (same conditions as for Fig. |5]l. The grid has 
= 1000 and M = 100 (signed ^-variable). 

sembling Eq. |2]provided the modulus k is changed for the "soft- 
ened modulus": 



2 



Vfl2 + + 2aR + 



(21) 



where A is called the "softening length" (generally considered as 
a free parameter). They found that the best expression for A that 
preserves the Newtonian character of the disk potential takes the 
form 



Ak hx g{a,R, . . .), 



(22) 



where h is the local semi-thickness of the disk, and g is a 
bounded, slowing varying function of the radius. This result is 
limited to disks with small aspect ratios, locally, i.e. to geomet- 
rically thin disks in the classical sense. This result has confirmed 
the common idea that /I is a certain fraction of the disk thickness. 
In the absence of mass density gradients in the direction per- 
pendicular to the disk equatorial plane, g goes through a sharp 
minimum of e"' k 0.368 ?AR - a (i.e. for u - 1). The sensitiv- 
ity with stratification appears weak since for a vertical parabolic 
profile, this minimum is e ~ 0.264 and the sharpness of the 
peak is unchanged. Appendix ICl reproduces the result prelimi- 
narily derived in the limit of g far from kernel singularity (i.e. for 
M — > or oo). In particular, we find that g tends asymptotically to 
3"'^^ ~ 0.577 in the homogeneous case, while it is 5"'''^ ~ 0.447 
in the parabolic case. We conclude that, for a given vertical pro- 
file, the function g does not vary much, except in the vicinity of 
R - a where it decreases by a factor of ~ 37%. As A does not 
depend on the pr ecise s hape h{a) provided the disk has a small 
aspect ratio (see iHure & Pierens 2009), this conclusion seems 
quite robust. This point is confirmed by the numerical experi- 
ments. In the following, we therefore neglect any radial variation 
of g, and consider that it is a constant whose nominal value will 
be selected a posteriori (and this assumption will be justified). 
Now, reversing Eq. |2T]as done before and using Eq. |22] we find 
at a given colatitude 



2m 



+ 1 + cotan^e - 0, 



(23) 
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where we have set e - h/a. It means that, provided the aspect 
ratio of the disk e is a constant, the ODE can still be formed. 
After some algebra and with the same non-dimensionalization 
as considered above, we actually find 



-— - (1 + s)— -5i = 0, 

ar r 



where: 



1 



[ V Mout^.i,outK(A;^_out) 



(24) 



(25) 



where = A:^ at a = Oin, and A;^_out = A;^ at a = aout- We see that 
S A only differs from Eq.[T4]by the presence of softened modulus 
A:^ (instead of k). Since we have 



k, < 



1 



+ e 



2„2 



< 1, 



(26) 



it follows that 5^ is no longer singular at the disk edges. As a 
result, the radial acceleration at the edges becomes finite. The 
role of A is to mimic an extra dimension (the vertical one). 

The softened potential t}/,\ entering Eq.|24]can then be solved 
numerically in the whole space, as done for the flat disk. In this 
process, the inner boundary value is slightly lower and must be 
changed to 



1 



< 1. 



(27) 



By looking at Eqs.[8]and|24l we see that there is a perfect 
continuity between these two ODEs (with and without softening 
length) as /I — » 0. It means that ^ i/r as /I — > (including 
the value at the origin). Also, it is worth noting that the explicit 
dependence with A must be regarded as a dependence with the 
disk thickness. 

To check the reliability of this approach, we first compared 
the softened potential i^ j determined numerically from Eq. |24] 
with the Newtonian potential of a geometrically thin disk with 
the same edges, same surface density (and mass), and semi- 
thickness h cc a (precisely e - 0.1) (hereafter thin disk configu- 
ration A^^This reference is computed from the splitting method 
(lHurell2005l) . which is very accurate. A typical result, limited to 
the vicinity of the inner edge, is shown in Fig. [T] We see that 
the agreement between the two curves is very good. There is a 
very weak sensitivity to the g-parameter. The best agreement is 
obtained for g ^ 0.486, which is very close to the average of 
bounds. Interestingly enough, thickness effects mainly shift the 
potential curve, leaving the radial gradients almost unchanged 
(except at the edges). Figure|8]shows the logarithm of the relative 
error in the entire midplane. The error is low, less than 0.1% in- 
side the disk. With g - -q^i the agreement is even better at a large 
distance from the disk (but not inside it). Figure |9] displays the 
error for all radii and all colatitudes. We conclude that the global 
solution of the softened ODE is very close to the Newtonian po- 
tential of a thin disk. The largest deviations are observed in the 
neighborhood of the disk surface (i.e. Z « h). This is expected 
since the softened modulus has been approximated here (but this 
point can be revisited if necessary). 

Up to now, only disks with a constant aspect ratio and a ho- 
mogeneous vertical mass density profile have been considered. 
The generality of the method can be studied easily by consid- 
ering other disk configurations, in particular, configurations for 



I 0.9 
o 



0.8 



0.486 




reference values (splitting method) 
ODE (flat disc, i.e. X=0) 
"softened" ODE 



0.1 



0.2 



0.3 



0.4 



0.5 



Fig. 7. Comparison of different potentials: Newtonian potential 
of a geometrically disk with e = 0.1 {red circles), Newtonian 
potential of a flat disk {dashed line), and softened potential with 
g - {i, 0.486, {thin lines). Values are normalized to their 
central value. In all cases, A = 0.1, Oout/flo = 1, and s - -1.5 
(and the same disk mass). 




Fig. 8. Relative deviation between the softened potential deter- 
mined in the midplane from Eq. |24]with g - 0.486, and 

the Newtonian potential of a geometrically thin disk with same 
edges and same mass (see also Fig.|7]i. 



which the conditions for the derivation of the ODE are violated. 
We illustrate this point by considering two other thin disk con- 
figurations: 

- configuration B: same as configuration A but the disk is 
strongly flared, with hi a - e(a/aout)°'^^ (this is an extreme 
case if we refer to most disk models). In this case, the expres- 
sion for the softening length is still valid, but the conditions 
for deriving the softened ODE are violated since A/ his not a 
constant. 
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Fig. 9. En-or map. Same as for Fig. |8]but at all colatitudes. 



- configuration C: same as configuration A but the mass den- 
sity profile is parabolic, vertically. In this case, the ODE is 
fully valid, but g must be set to the appropriate value (in be- 
tween e"'*/^ and 5"'^^; see above). 

In all cases, the disk is geometrically thin at the outer edge (with 
e - 0.1), which is required to keep the softening length valid). 
We performed the same comparisons as for configuration A and 
the results are summarized in Fig. [10] We find that the relative 
deviation between the Newtonian potential of the thin disk and 
the softened potential (found from the softened ODE) typically 
varies as shown in Fig. [8] and is much less than 1% inside the 
disk. This is also the case outside the disk, except for configura- 
tion B for which the relative error gradually rises from R = Oout 
and reaches 10% at infinity. We therefore conclude that using 
the softened ODE to mimic the Newtonial potential of a thin 
disk with a precision of a few percent is fully justified as soon 
as the disk is geometrically thin at all radii. The best values are 
reproduced when the flaring angle of the thin disk is close to 
zero. 



6. Concluding remarks 

In this article, we have shown that the gravitational potential of a 
flat, power-law disk can be numerically determined in the whole 
physical space by solving, for any colatitude, a linear system 
resulting from the second-order discretization of a first-order or- 
dinary diff'erential equation (ODE). The computing time is es- 
pecially low because it is proportional to the number of mesh 
points. It is also demonstrated that the softened potential based 
on the prescription for the softening length by iHure & PierensI 
(120091) obeys a similar ODE whose solution agrees very well 
with the potential of a geometrically thin disk with same size, 
mass, and edges and a constant aspect ratio. We briefly tackle 
the limits of the method, and conclude that the Newtonian po- 
tential is reproduced in the full space, provided i) the disk is 
geometri cally thin at all radii, which justifies the use of A pro- 
posed bv lHure & PierensI (l2009l) . and ii) the local flaring remains 
moderate (i.e. h/a does not vary much with the radius), hence 
supporting the derivation of the ODE. The implementation of 
a softening length not only mimics a vertical stratification, but 



— - flat disk (i.e. X=0) 
o geometrioaiiy tliin disk 
— "softened" ODE 



strong flaring 




with X =]/3 



variable ? 

Fig. 10. Same legend as for Fig.|8]but for configuration A, B, and 
C. For clarity, potential values are shifted upward for configura- 
tion B, and downward for configuration C. 



also totally removes edge singularities typical of flat disks with 
sharp edges. If the surface density does not obey a power law 
but can be written in the form of a series of power laws (e.g. 
Taylor expansion), the potential is then obtained easily by super- 
position (see Paper I). The computing time is then proportional 
to the number of terms in the series. This approach is especially 
efficient at generating grids of forces. This offers reliable tools 
for investigating the dynamics of particles in a system contain- 
ing a perturbing, gravitating thin disk (e.g. ISubr & Karasll2005l) 
and possibly for linking back various families of trajectories to 
the mass distribution in the disk. The two components of the 
gravitational force are found from potential values by finite dif- 
ferences. The radial component can, however, be found directly 
from the ODE. For a unit mass, this component is simply 



F, = -(1 + s) 



A, 



(28) 



It would be interesting to investigate the existence of a first-order 
differential equation associated with the colatitudinal gradient of 
the potential, complementary to the radial ODE reported here. 
This is the subject of ongoing work. 
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Appendix A: Derivation of the generalized ODE 

Starting from Eq. |7] and assuming 9 constant, the derivative of u 
with respect to k is 

du 8m^ sin^ 6 



dk -M^sin^e) 



du 
dk 



(A.l) 



This derivative is therefore a function of only k. The potential 
in Eq. |2]can then be written as the product of a function of the 
spherical radius r by an integral over k, namely 



where 



r 



'Hik)dk, 



Hik) = m'+2K(A;)A:x — . 

dk 



(A.2) 



(A.3) 



The new integral bounds, k\a and A;out, correspond to the disk 
edges ain and Oout, and they are found from Eq. |4] In particular, 
we have 



Aa-,„R 



+ + 2a\nR 



(A.4) 



and 



k^ - 

"'nut 



(A.5) 



respectively. For constant colatitude 6, the derivative of tfr with 
respect to the spherical radius is then given by 



di// 
dr 



(1 + .9)^-2GIoV/?'' 



(A.6) 



'2t.\ 

dr 



This quantity is just the opposite of the radial acceleration due to 
the disk. Using an elementary derivation rule (see Paper I), we 
can calculate the right hand-side of this equation. We find 



dr r 



where 

dko^ii _ dk 



dr 



{A.l) 



(A.8) 



and similarly for ^in. When a is held constant, we have 



dk 
dr 

and so 
dk 
d? 



k^ [u^sin^e- l) 



8a sin 



du 
dk 



r^ sin 6 



It follows that 

?-(i + .)!^-A = o, 

dr r 
where A is defined by: 



r 



(A.9) 



(A. 10) 



(A. 11) 



(A. 12) 
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Mout - aout/R and - ain/R- This ordinary differential equa- 
tion (ODE) is the generalization to the whole space of the ODE 
reported in Hure & Hersant (2007), which was valid only in the 
disk plane (i.e. the case 6 = |). This expression differs mainly 
by the presence of the spherical radius r (instead of R) and by 
the presence of the Z-dependent modulus k (instead of m). As in 
Paper I, the second member A is an analytical function of a sin- 
gle spatial coordinate r. It now depends on the colatitude 6 and 
on the four parameters s, ao, ai„, and Oout- 



Appendix B: Asymptotic properties 

The differential equation possesses the right properties both at 
small and at large distances. At a large distance from the disk, 
the modulus in the complete elliptic integrals tends to zero. We 
have 



limK(;t) = -, 



(B.l) 



and then 



lim S 



(2 + s)x2^ 



Since the total mass of the disk is given by the expression 



(B.2) 



Md = 271 I 'L(a)ada = _ZllZ^{L±li^ (B.3) 
the above ODE can then be rearranged into 

^(l+.)^ + (2 + .)e^ 
,,(^+c^) + ^ + 2^ (B.4) 

As this equation must be satisfied for any s, we must have 
GMd 



lim tfr 



which is the expected behavior (the disk is no longer distinguish- 
able from a point mass). This also implies that 



# GMd 
lim -— = —. 

r->co dr r 



(B.6) 



At a short distance around the origin (i.e. r <K Uin), we can 
perform a second order expansio n of the S -term by expandin g 
the elliptic integral accordingly (iGradshtevn & RvzhikI Il965l) . 
We find 



lim rS - -(1 -1- i) - Mqu, . 



and then the ODE becomes 
ar 

whose solution is 

1 , -V^-l -2 



Aq sin^ 6 _2 



(B.7) 



(B.8) 



(B.9) 



At second order, the potential in the inner domain (r <K ain) is 
quadratic with the cylindrical radius R (while the gravitational 
acceleration is linear). 



Appendix C: The softening lengtli at large relative 
separations 

As in lHure &Pierensl(l2009h . we consider a vertical stratification 
of the form 



p(z) = po 



1 



(C.l) 



for |z| < h (and elsewhere), where po is the density at the disk 
midplane, h the local semi-thickness (both a function of the ra- 
dius a in general), and ^ > 1 an integer Homogeneous profiles 
correspond to ^ oo, whereas the parabolic profile is obtained 
for q = 1. We can compute the contribution of vertical stratifica- 
tion to the potential from the integral 



where 

Xcj = i2q + 1) 



X'(i) 



2q J 

kK(k)d-. 
h 



(C.2) 



(C.3) 



The approxima t ions k x 1 and K{k) ^ In4/A:', assumed in 
[Hure & Pierens' (2009), is only valid for a ^ R. To derive the 
asymptotic behavior at large relative separation, it is suflicient to 
consider the case — > 0, that is, at second-order 



-1—1 

2\a+R) 



as well as K(;t) * f . We then find 



and so 



(B.5) X 



2q + 1 
2(2^ + 3) \a+R 



2q + 1 



6(2^ -H 3) \a+R 



(C.4) 



(C.5) 



(C.6) 



The softening length is then deduced by equating - ms|, 
where is the softened modulus: 



2 Vfl^ 



V(fl + R)- + 
We then have 



2g+l 
'2(2^-h3)' 



(C.7) 



(C.8) 



